Fragmentation processes in impact of spheres 
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We study the brittle fragmentation of spheres by using a three-dimensional Discrete Element 
Model. Large scale computer simulations are performed with a model that consists of agglomerates 
of many particles, interconnected by beam-truss elements. We focus on the detailed development of 
the fragmentation process and study several fragmentation mechanisms. The evolution of meridional 
cracks is studied in detail. These cracks are found to initiate in the inside of the specimen with 
quasi-periodic angular distribution. The fragments that are formed when these cracks penetrate 
the specimen surface give a broad peak in the fragment mass distribution for large fragments that 
can be fitted by a two-parameter Weibull distribution. This mechanism can only be observed in 
3D models or experiments. The results prove to be independent of the degree of disorder in the 
model. Our results significantly improve the understanding of the fragmentation process for impact 
fracture since besides reproducing the experimental observations of fragment shapes, impact energy 
dependence and mass distribution, we also have full access to the failure conditions and evolution. 

PACS numbers: 46.50.-fa, 62.20.Mk, 05.10.-a 



I. INTRODUCTION 

Comminution is a very important step in many indus- 
trial applications, for which one desires to reduce the 
energy necessary to achieve a given size reduction and 
minimizing the amount of fine powder resulting from 
the fragmentation process. Therefore a large amount 
of research has already been carried out to predict the 
outcome of fragmentation processes. Today the mecha- 
nisms involved in the initiation and propagation of single 
cracks are fairly well understood, and statistical models 
have been applied to describe macroscopic fragmenta- 
tion [TJ [5]. However, when it comes to complex frag- 
mentation processes with dynamic growth of many com- 
peting cracks in three-dimensional space (3D), much less 
is understood. Today computers allow for 3D simula- 
tions with many thousands of particles and interaction 
forces that are more realistic than simple central poten- 
tials. These give a good refined insight of what is really 
happening inside the system, and how the predicted out- 
come of the fragmentation process depends on the system 
properties. 

Experimental and numerical studies of the fragmenta- 
tion of single brittle spheres have been largely applied to 
understand the elementary processes that govern com- 
minution laiiigEiiziiHiEiEiiiiiiiiiiaiiiiiiiiiniiizi 

HSllIllEOlEIlES]- Experiments that were carried out in 
the 60s analyzed the fragment mass and size distributions 
[31 m [S] with the striking result that the mass distribu- 
tion in the range of small fragments follows a power law 
with exponents that are universal with respect to mate- 
rial, or the way energy is imparted to the system. Later 
it became clear that the exponents depend on the di- 
mensionality of the object. These results were confirmed 
by numerical simulations that were mainly based on Dis- 
crete Element Models (DEM) [151 HI [551 HHj- For 



large fragment masses, deviation from the power law dis- 
tribution could be modeled by introducing an exponential 
cut-off, and by using a bi-linear or Weibull distribution 
[El [m [m Ea EH EH So]. Another important finding 
was, that fragmentation is only obtained above a certain 
material dependent energy input [51 [SI [31]. Numerical 
simulation could show that a phase transition at a crit- 
ical energy exists, with the fragmentation regime above, 
and the fracture or damaged regime below the critical 
point [IS|[Tni[lT]. 

The fragmentation process itself became experimen- 
tally accessible with the availability of high speed cam- 
eras, giving a clear picture on the formation of the frag- 
ments |5l[6l[Il[8l[9l[10l[IIl[T2l[13]. Below the critical 
point, only slight damage can be observed, but the spec- 
imen mainly keeps its integrity. Above but close to the 
critical point, the specimen breaks into a small number of 
fragments of the shape of wedges, formed by meridional 
fracture planes, and additional cone-shaped fragments at 
the specimen-target contact point. Way above the criti- 
cal point, additional oblique fracture planes develop, that 
further reduce the size of the wedge shaped fragments. 

Numerical simulations can recover some of these find- 
ings, but while two-dimensional simulations cannot re- 
produce the meridional fracture planes that are respon- 
sible for the large fragments [H [SI HIl US EHl EI] , three- 
dimensional simulations have been restricted to relatively 
small systems, and have not focused their attention on 
the mechanisms that initiate and drive these meridional 
fracture planes |15L 119] . Therefore, their formation and 
propagation is still not clarified, although the resulting 
two to four spherical wedged-shaped fragments are ob- 
served for a variety of materials and impact conditions 
in [3 [m [20] . Arbiter et al.| |5j argued, based on the anal- 
ysis of high speed photographs, that fracture starts from 
the periphery of the contact disc between the specimen 
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and the plane, due to the circumferential tension induced 
by a highly compressed cone driven into the specimen. 
However, their experiments did not allow access to the 
damage developed inside the specimen during impact. 
Using transparent acrylic resin, ,Majzoub and Chaudhri| 
[8] observed damage initiation at the border of the con- 
tact disc, but in their experiments plastic flow and ma- 
terial imperfections may have a dominant role. 

In this paper we present three-dimensional simulations 
of brittle solid spheres under impact with a hard plate. 
With our simulations, the time evolution of the fragmen- 
tation process and stress fields involved are directly ac- 
cessible. We have focused our attention on the processes 
involved in the initiation and development of fracture, 
and how they lead to different regimes in the resulting 
fragment mass distributions. Our results can reproduce 
experimental observations on fragment shapes, impact 
energy dependence, and mass distributions, significantly 
improving our understanding of the fragmentation pro- 
cess in impact fracture. 



II. MODEL AND SIMULATION 
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of the contact force include damping and friction forces 
and torques in the same way as described in Refs. 

[HIIHllllIlD]. 



r.. 













Discrete Element Models (DEM) have been success- 
fully u sed since they were introduced by [Cundall and 
Strack to study rock mechanics [32]. Applications range 
from static, to impact and explosive loading, using 
elementary particles of various shapes that are con- 
nected by different types of massless cohesive elements 
[llinillHlIiniEOliaESllMllSgE^ in gen- 

eral, Newton's equation governs the translational and 
rotational motion of the elements, that concentrate the 
whole mass. Forces and torques arise from element in- 
teractions, from the cohesive elements, volumetric forces, 
and of course from interaction with boundaries like walls. 

Throughout this work we use a three-dimensional (3D) 
implementation of DEM where the solid is represented 
by an assembly of spheres of two different sizes. They 
are connected via beam-truss elements that deform by 
elongation, shear, bending, and torsion. The total force 
and moment acting on each element consists of the 
contact forces resulting from sphere-sphere interactions, 
pc _ pov _|_ pdiss ^ ^Yic stretching and bending forces 

= F'^^° + Q and moments transmitted by the 
beams attached. 

The contact force has a repulsive term due to elastic in- 
teraction between overlapping spherical elements, which 
is given by the Hertz theory [35] as a function of the ma- 
terial Young's modulus , the Poisson ratio , and the 
deformation ^. The force on element j at a distance r^j 
relative to element i (see Fig. [ija) ) is given by 
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where the overlapping distance 



Ri 



describes the deformation of the spheres, 
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FIG. 1: (a) Representation of the overlap interaction between 
two elements, (b) Typical deformation of a beam in the x- 
y plane, showing the resulting bending and shear forces and 
torques. The z-axis is perpendicular to the image. 



The 3D representation of beams used in this work is an 
extension of the two-dimensional case of Euler-BernouUi 
beams described in Ref. [4lj. In 3D, however, the total 
deformation of a beam is calculated by the superposition 
of elongation, torsion, as well as bending and shearing 
in two different planes. The restoring force acting on 
element j connected by a beam to element i due to the 
elongation of the beam is given by 



(2) 



where E^ is the beam stiffness, £ = {\rij \ — lo) /Iq, with 
the initial length of the beam and its cross section . 

The flexural forces and moments transmitted by a 
beam are calculated from the change in the orientations 
on each beam end, relative to the body-fixed coordinate 
system of the beam (e^,ey,e^). Figure mb) shows a ty- 
pical deformation due to rotation of botn beam ends rel- 
ative to the axis, with e\ oriented in the direction of 
rtj. Given the angular orientations Of and the cor- 
responding bending force Q^'^ and moment M^'^ for the 
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elastic deformation of the beam are given by [IT] : 




where / is the beam moment of inertia. Corresponding 
equations are written for general rotations around Cy , and 
the forces and moments are added up. Additional torsion 
moments are added to consider a relative rotation of the 
elements around e^: 

M-'b ^ -GH'°'-^-^—^ei, (4) 
■' L 

with and 7*°'" representing the shear modulus and mo- 
ment of inertia of the beams along the beam axis, respec- 
tively. The bending forces and moments are transformed 
to the global coordinate system before they are added to 
the contact, volume and walls forces. 

Beams can break in order to explicitly model damage, 
fracture, and failure of the solid. The imposed breaking 
rule takes into account breaking due to stretching and 
bending of a beam |M1 SOI SI] , which breaks if 

/ £ y ^ max(|0,|,|0,|) 

\£thj Qth ~ 

where e = AZ/Zq is the longitudinal strain, and 9^ and 
0j are the general rotation angles at the beam ends be- 
tween elements « and j , respectively. Here cos di — e!^'' • , 
where {e'tj^ , e^y , ef) define the i— particle's orientation in 
the beam body-fixed coordinate system, similar calcula- 
tion is performed to evaluate 9j. Equation ^ has the 
form of the von Mises yield criterion for metal plasticity 
[40] |43] . The first part of Eq. ^ refers to the breaking 
of the beam through stretching and the second through 
bending, with eth and 9th being the respective thresh- 
old values. The introduced threshold values are taken 
randomly for each beam, according to the WeibuU distri- 
butions: 



P{eth) = 
P{Oth) = 




Here k , Eq and 9^ are parameters of the model, control- 
ling the width of the distributions and the average values 
for Eth and 9th respectively. Low disorder is obtained by 
using large k values, large disorder by small k. Disor- 
der is also introduced in the model by the different beam 
lengths in the discretization as described below. 

The time evolution of the system is followed by nu- 
merically solving the equations of motion for the transla- 
tion and rotation of all elements using a 6"'-order Gear 
predictor-corrector algorithm, and the dynamics of the 



rotations of the elements is described using quaternions 
[41] |44] . The breaking rules are evaluated at each time 
step. The beam breaking is irreversible, which means 
that broken beams are excluded from the force calcula- 
tions for all consecutive time steps. 



System formation and characterization 

Special attention needs to be given to the discretiza- 
tion in order to prevent artifacts arising from the sys- 
tem topology, like anisotropic properties, leading to non 
uniform propagation of elastic waves or preferred crack 
paths. In our procedure we first start with 27000 spher- 
ical elements that we initially place on a cubic lattice 
with random velocities. The element diameters are of 
two different sizes, with D2 = 0.95Di, that are randomly 
assigned, leading to more or less equal fractions. Once 
the elements are placed, the system is left to evolve for 
50000 time steps, using periodic boundary conditions, in 
a volume that is about 8 times larger then the total vol- 
ume of the elements. This way we obtain truly random 
and uniformly distributed positions. 

To compact the elements, a centripetal constant ac- 
celeration field, directed towards the center of the sim- 
ulation box, is imposed. Due to this field the elements 
form a nearly spherical agglomerate at the center of the 
box. The system is allowed to evolve until all particle 
velocities are reduced to nearly zero due to dissipative 
forces. 

With the elements compacted, the next stage is to con- 
nect them by beam-truss elements. This is achieved in 
our model through a Delaunay triangulation of their po- 
sitions. As a consequence, not only spherical elements 
that are initially in contact or nearly in contact with 
each other are connected, but the resulting beam lat- 
tice is equivalent to a discretization of the material us- 
ing a dual Voronoi tessellation of the material domain 
[43] |45j |46] . After the bonds have been positioned, their 
Young's moduli are slowly increased while the centripetal 
field is reduced to zero. During this process the material 
expands to an equilibrium state, reducing the contact 
forces. The bond lengths and orientations are then reset 
so that no initial residual stresses are present in the beam 
lattice. The final solid fraction obtained is approximately 
0.65. We have compared impact simulations of specimens 
compacted as described above with specimens using ran- 
dom packings of spheres as reported in Ref. [17] , which 
have no preferential direction in the packing process such 
as the one that could be imposed by the centripetal field. 
No significant difference was found in the simulation re- 
sults, indicating that possible radially aligned locked-in 
force chains are not relevant. 

Once the system is formed, the specimen is shaped to 
the desired geometry by removing particles and beams 
that are situated outside the chosen volume. The mi- 
croscopic properties, namely the elastic properties of 
the elements and bonds, as well as the bond breaking 
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thresholds, are chosen to attain the desired macroscopic 
Young's modulus, Poisson's ratio, as well as the tensile 
and compressive strength. Table 1 summarizes the in- 
put values used in the simulations presented in this pa- 
per. These were chosen to obtain macroscopic proper- 
ties close to the mechanical properties of polymers like 
PMMA, PA, and nylon at low temperatures. Figure [2] 
displays the stress-strain curve measured by quasi-static, 
uni-axial tensile loading of a bar, as depicted in the inset. 
The microscopic and resulting macroscopic properties are 
resumed in Table 1, for a sample size (16 x 8 x 8 mm). 
The experiment is performed by measuring the force re- 
quired to slowly move the upper and lower surfaces (see 
inset of Fig.[2| at a constant strain rate of 0.004 s^^. The 
stress-strain curve is basically linear until the strength is 
reached where rapid brittle fracture of the material takes 
place. Oscillations in the broken specimen fractions can 
be seen after the system is completely unloaded due to 
elastic waves. The Young's modulus measured from the 
slope of the curve is 7.4±0.5GPa, is presented along with 
other macroscopic properties of the material in Table 1. 



in the negative z-direction, assigned to all its composing 
elements. The computation continues until no additional 
bonds are broken for at least 50 us. 

For comparative reasons we calculate the evolution of 
the stress field using an explicit Finite Element (FE) 
analysis. The FE model is composed of axisymmetric, 
linear 4-node elements with macroscopic properties taken 
from the results of the DEM simulations (see Table 1). 
Along the central axis through the sphere and ground 
plate, symmetry boundary conditions are imposed, the 
bottom of the target plate is encastred and contact sur- 
faces for the sphere and plate are defined. Figure [3](a) 
shows a comparison between the impact simulation using 
our DEM model and a Finite Element Model simulation. 
In Fig.jsja), the DEM elements are colored according to 
the amplitude of their accelerations to show the propaga- 
tion of a longitudinal shock wave that was initiated at the 
contact point. The wave speed can be estimated to be ap- 
proximately 2200 ± lOOm/s, which is consistent with the 
Young's modulus of the material derived from Fig. [2] and 
its density. The time evolution of the potential energy 
stored in the system is compared in Fig. [3](b), showing 
excellent quantitative agreement between the two mod- 
els. 

After the characterization of the system properties we 
allow for the cohesive elements to fail in order to study 
the fragmentation properties. 



n.if) 
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FIG. 2: Stress-strain curve for specimen under quasi-static 
loading. The inset shows the load geometry. Abrupt brittle 
fracture behavior can be observed at about e = 0.019. 

In order to simulate the impact of a sphere on a fric- 
tionless hard plate, a spherical specimen with diame- 
ter D — 16 mm is constructed, and a fixed plane with 
Young's modulus 70 GPa is added to the simulation. The 
spherical specimen has a total of approximately 22000 el- 
ements, with around 32 across the sample diameter. The 
contact interaction between the elements and the plate 
is identical to the element-element contact interaction, 
only with £, ^ Ri — rip , where r^p is the distance be- 
tween the particle center and the plate. The specimen 
is placed close to the plate with an impact velocity Vi, 




time (,,,!) 



FIG. 3: (a) Comparison of deformations and shock-wave 
propagation obtained between DEM and FEM simulations 
for Vi — 117 m/s. (b) Time evolution of the elastic potential 
energy stored in the system for the same velocity obtained by 
DEM (solid line) and FEM (dashed line) simulations. 
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TABLE L Micro- and macroscopic material model properties. 
Typical model properties (DEM): 



Beams: 



stiffness 


E*'/G*' 


6 


GPa 


average length 


L 


0.5 


mm 


diameter 


d 


0.5 


mm 


strain threshold 


eo 


0.02 


_ 


bending threshold 




3 


o 


shape parameter 


k 


3/10 


- 


Particles: 








stiffness 




3 


GPa 


diameter 


Di 


0.5 


mm 


density 


P 


3000 


kg/m^ 


Hard plate: 








stiffness 




70 


GPa 


Interaction: 








friction coefficient 


M 


1 


- 


damping coefficient 




0.25 


s-i 


friction coefficient 


It 


0.05 


s"^ 


System: 








time increment 


At 


le-8 


s 


number of particles 


NP 


22013 




number of beams 


N'' 


135948 


- 


solid fraction 




0.65 




sphere diameter 


D 


16 


mm 


Macroscopic properties 


1 (DEM): 




system stiffness 


E 


7.4 ±0.5 


VjrJT cL 


Poisson's ratio 


V 


0.2 




density 


P 


1920 


kg/m^ 


system strength 




110 


MPa 


Comparison: 








DEM 




FEM 




longitudinal 2210 ±100 


2270 ± 20 


m/s 


wave speed 








contact time 31.4 




31.4 


^s 



III. FRAGMENTATION MECHANISMS 

In this section we explore the different fragmentation 
mechanisms in the order of occurrence and increasing 
energy input. The first yield that arises in the material 
is diffuse damage that occurs in the region above the 
contact disc. It can be seen from Fig. |4]ja), that this 
damage region is centered in the load axis, at a distance 
approximately I?/4 from the plane. 

We can see a strong correlation of the position of the 
diffuse cracking in the DEM results (Fig. |4]^a)) with the 
location of a region with a bi-axial stress state in the x-y- 
plane and a superimposed compression in the z-direction, 
as calculated using FEM (Fig. gb)), also in agreement 
with experimental results reported in Ref. [B] . This re- 
sult, along with the one presented in Fig. [3] suggests that 
the use of three-dimensional beams, as compared with 
the use of simple springs, despite of the reduced number 
of degrees of freedom in the breaking criterion, could re- 



cover quite well the influence of complex stress states in 
the crack formation in a more precise way. 




FIG. 4: Initial damage due to bi-axial stress state, (a) Ver- 
tical cut through the center of the sphere from the DEM sim- 
ulation showing broken bonds represented by dark color, (b) 
Stress fields calculated with FEM model. Left side are shear 
stresses in global coordinates from to -400MPa (black to 
white) while the right side shows circumferential stresses in 
local spherical coordinates with the center in the sphere center 
ranging from to 130MPa (black to white). 

As time evolves, meridional cracks start to appear. 
The origination of this type of cracks is explored in Fig.js] 
where we plot in Fig. [sj^a) the positions of the broken 
bonds in two different projections, showing well defined 
meridional crack planes that propagate towards the lat- 
eral and upper free surfaces of the specimen. In Fig. |5]^b) 
we plot the angular distribution of the broken bonds for 
different times. Here g{ff) is the probability of finding 
two broken bonds with an angular separation Q. Note 
that their positions are projected into the plane perpen- 
dicular to the load axis. The evident peaks in g{Q) are 
a clear indication that the cracks are meridional planes 
that include the load axis. In this particular case, the 
cracks are separated by an average angle of about 60 de- 
grees, and they become evident 13 to 15 /is after impact 
(v, = 120 m/s). 

In order to understand what governs the orientation 
and angular separation of these meridional cracks we per- 
formed many different realizations with different seeds of 
the random number generator and impact points. For 
all cases the orientation of the cracks can change, but 
not their average angular separation. We observe that 
for strong disorder (Eq. (|6|), a larger amount of uncorre- 
lated damage occurs, but the average angular separation 
of the primary cracks does not change. This suggests 
that the formation of these cracks arises due to a com- 
bination of the existence of local disorder and the stress 
field in the material, but does not depend on the degree 
of disorder. 

As we can see from the FE calculations and from the 
damage orientation correlation plot (Fig. [s]) inside of the 
uniform biaxial tensile zone, no preferred crack orienta- 
tion is evidenced. Many microcracks weaken this ma- 
terial zone, decreasing the effective stiffness of the core. 
Around the weakened core the material is intact and un- 
der high circumferential stresses. It is in this ring shaped 
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FIG. 5: (a) Colored dots display the positions of the broken 
bonds according to the time of breaking, (b) Angular distri- 
bution function of broken bonds as a function of the angular 
separation when their positions are projected into the plane 
perpendicular to the load axis. 



zone, that we observe to be the onset of the meridional 
cracks when we trace them back. For increasing impact 
velocity we observe a decrease in the angular separation 
of crack planes and thus more wedge-shaped fragments. 
Therefore this fragment formation mechanism can not be 
explained by a quasi-static stress analysis. The observa- 
tion is in agreement with experimental findings and can 
be explained by the basic ideas of Mott's fragmentation 
theory for expanding rings [48 . Due to the stress release 
front for circumferential stresses, once a meridional crack 
forms, stress is released in its neighborhood; the fractured 
regions spread with a constant velocity and the probabil- 
ity for fracture in neighboring regions decreases. On the 
other hand in the unstressed regions, the strain still in- 
creases, and so does the fracture probability along with it. 
The average size of the wedge shaped fragments therefore 
is determined by the relationship between the velocity of 
the stress release wave and the rate at which cracks nucle- 
ate. Thus the higher the strain rate, the higher the crack 
nucleation rate and the more fragments are formed. We 
measured the strain rate at different positions inside the 
bi-axially loaded zone, finding a clear correlation between 
impact velocity and strain rate. Even though we frag- 
ment a compact sphere and not a ring, when it comes to 
the formation of meridional cracks, we observe that they 
form in a ring shaped region and that Mott's theory can 
qualitatively explain the decrease of angular separation 



between wedge shaped fragments with increasing impact 
velocity. 

If enough energy is still available, some of the merid- 
ional plane cracks grow outwards and upwards, breaking 
the sample into wedge shaped fragments that resemble 
orange slices. 

As the sphere continues moving towards the plate, a 
ring of broken bonds forms at the border of the con- 
tact disc due to shear failure (compare Fig.[6][a)). When 
the sphere begins to detach from the plate, the cone has 
been formed by high shear stresses in the contact zone 
(see Fig. ^h) left) by a ring crack that was able to grow 
from the surface to the inside of the material under ap- 
proximately 45° (Fig. |6[b)). It detaches, leaving a small 
number of cone shaped fragments that have a smaller 
rebound velocity than the rest of the fragments due to 
dissipated elastic energy (Fig. [6jc)). 



— fmWWWimM immMmkmlm immmmwmm 



FIG. 6: Vertical meridional cut of the sphere at different 
stages during impact, showing the separation of lower frag- 
ments, (a) Formation of a ring of broken bonds due to shear 
failure, (b) These broken bonds evolve into cracks that prop- 
agate inside the sample, (c) Finally these cracks lead to the 
detachment of the lower fragments. 

Oblique plane cracks may still break the large frag- 
ments further, if the initial energy given to the system is 
high enough. Therefore they are called secondary cracks. 
Figure |7|a) shows a vertical meridional cut of a sample 
where these cracks can be seen. The intact bonds are 
colored according the final fragment they belong to. 




FIG. 7: (a) DEM simulation at Vi = 140 m/s exemplifying 
the secondary cracks. The bonds are colored according to 
the final fragment they belong to. (b) 2D simulations using 
polygons as elementary particles [3T]. 
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These secondary cracks are very similar to the oblique 
cracks observed in 2D simulations yL4J 121] . For compari- 
son we show in FigjTj^b) the crack patterns obtained from 
a 2D DEM simulation that uses polygons as elementary 
particles. In the 2D case, we observe a cone of numerous 
single element fragments and meridional cracks cannot 
be observed. 



IV. FRAGMENTATION REGIMES 



The amount of energy necessary to fragment a material 
is a parameter that is very important for practical appli- 
cations in comminution. In fragmentation experiments 
two distinct regimes for damage and fragmentation can 
be identified depending on the impact energy: below a 
critical energy [H |6l [19] damage takes place, while above 
fragments are formed. Figure |8] compares the final crack 
patterns after impact with different initial velocities. The 
intact bounds are colored according to the final fragment 
they belong to, and gray dots display the positions of bro- 
ken bonds. The fragments have been reassembled to their 
initial positions to provide a clearer picture of the result- 
ing crack patterns. For the smaller impact velocities it is 
possible to identify meridional cracks that reach the sam- 
ple surface above the contact point, but fragmentation is 
not complete and one large piece remains (Figs, [sja) and 
[8][b)). We call these meridional cracks primary cracks, 
since as one increases the initial energy given to the sys- 
tem, some of them are the first to reach the top free 
surface of the sphere, fragmenting the material into a 
few large pieces, typically two or three fragments with 
wedge shapes (Fig. [Sj^c)). When we increase the initial 
energy secondary oblique plane cracks break the orange 
slice shaped fragments further (Fig. [sjd)). Additional 
increase in the impact velocity causes more secondary 
cracks and consequently the reduction of the fragment 
sizes (Figs.[8];e) and[8];f)). 

The shape and number of large fragments resulting 
from the numerical model for smaller impact energies, as 
well as the location and orientation of oblique secondary 
cracks for larger energies, are in agreement with experi- 
mental findings [TTlfniPU]. 

We can identify that for velocities smaller then a 
threshold value, the sample is damaged by the impact 
but not fragmented. This threshold velocity for frag- 
mentation has been found experimentally and numeri- 
cally [SJ H^l m]. In particular, it has been found from 
2D simulations that a continuous phase transition from 
damaged to fragmented outcome of impact fragmenta- 
tion can be tuned by varying the initial energy imparted 
to the system [15] [^. 

Following the analysis in references [TSJ |2T] the final 
state of the system after impact is analyzed by observing 
the evolution of the mass of the two largest fragments, 
as well as the average fragment size (shown in Fig. |9]). 
The average mass M2/M1, with = - M^^^ 

excludes the largest fragment. It can be observed that 




FIG. 8: Front view of the reconstructed spheres, showing 
the final crack patterns at the surfaces for different initial 
velocities. Gray dots are placed at the positions of broken 
beams while different colors are chosen for different fragments. 



below the threshold value vth — 115 m/s the largest frag- 
ment has almost the total mass of the system, while the 
second largest is orders of magnitude smaller. This be- 
haviour implies that for Vi < vtu the system does not 
fragment, it only gets damaged. For velocities larger then 
Vth the mass of the largest fragment decreases rapidly. 
The second largest and average fragment masses increase, 
having their maximum at 117.5 m/s for this material 
strength. 

The results shown in Fig. [9| are in very good qualita- 
tive agreement with those obtained from simulations in 
different geometries and load conditions [18] [STJ |49] , in- 
dicating that our model shows a phase transition from a 
damaged to a fragmented state. 
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that these fragments have their origin in mechanisms dis- 
tinct from the ones that form small fragments. Fig. 10 ^b) 




40 160 ISO 200 



FIG. 9: The mass of first and second largest fragment and the 
average fragment mass as a function of the impact velocity. 



RESULTING FRAGMENT MASS 
DISTRIBUTION 



One of the first and still most important character- 
izations for fragmentation processes are fragment mass 
distributions. Experimental and numerical studies on 
fragmentation show that the mass distribution follow a 
power law in the range of small fragments, whose expo- 
nent depends on the fragmentation mechanisms, while 
the mass distribution for large fragments is usually rep- 
resented by an exponential cut-off of the power law. The 
fragment mass distributions that are obtained from our 
three-dimensional simulations are given in Fig. 10 'a) for 
different impact velocities Vi. Here F{m) represents the 
probability density of finding a fragment with mass m 
between m and m + Am. Where m is the fragment 
mass normalized by the total mass of the sphere Mtot- 
The values are averaged over 36 simulations, changing 
the random breaking thresholds and randomly rotating 
the sample to obtain different impact points. For ve- 
locities below the critical velocity vtt of our model, the 
fragment mass distribution shows a peak at low frag- 
ment masses, corresponding to some small fragments. 
The pronounced isolated peaks near the total mass of 
the system correspond to large damaged, but still unbro- 
ken system (see also Figs. [Sja) and (b)). Fragments at 
intermediate mass range are not found for small initial 
velocities. At and above Vth, the fragment mass distri- 
bution exhibits a power law dependence for intermediate 
masses, F{m) ~ , (dashed line in Fig. 10 'a)) with 
T = 1.9 ± 0.2 inni EI], and a broad maximum can be 
observed in the histogram for large fragments, indicating 



shows the cumulative size distribution of the fragments 
weighted by mass, Q3, for the same velocities. Q3 is cal- 
culated by summing the mass of all the fragments smaller 
then a given size s. The size of a fragment is estimated as 
the diameter of a sphere with identical mass, the values 
are normalized by the sample diameter. By this repre- 
sentation the large fragments are better resolved. We 
can see that the shape of the size distribution for large 
fragments can be described by a two-parameter WeibuU 
distribution, Q3{s) = 1 — exp[— (s/sc)'^°] (dashed line in 
Fig. [I0];b), with = 0.75 and ks = 5.8). The WeibuU 
distribution is used here since it has been empirically 
found to describe a large number of fracture experiments, 
specially for brittle materials |29j . With increasing initial 
velocity, the average fragment size shifts towards smaller 
values, also in agreement with experimental findings from 

Refs. [niiaol. 



The local maximum in the fragment mass distribution 
for large fragments represents those fragments, that are 
formed by the meridional cracks. As we can observe from 
Fig. |ll[ the fragment mass distribution is independent of 
the amount of disorder or material that the specimen is 
composed of (fc is in the breaking thresholds distributions 
in Eq. Q). Near the critical velocity Vth we can iden- 
tify two main parts in the fragment mass distribution. 
For m up to around 1/40 (approximately 550 elements), 
the power law F{m) ^ mr"^ J{m/fho) with the cutoff 
function f{m/rho) containing an exponential component 
exp (—m/rho) can be used like in Ref. [49 . However, for 
larger m, F{m) can also be described by a two-parameter 
WeibuU distribution 



F{m) 



rhi 



m 

m 



exp 



m 
fhi 



(7) 



The dashed line in Fig. 11 corresponds to a fit to the 
data using TOo = 0.004 ±0.001, m = 0.3 ±0.02 and fc/ = 
1.9 ± 0.1. The good quality of the fit allows for a better 
estimation of the exponent of the power-law distribution 
in the small fragment mass range r — 2.2 ± 0.02. 

For the material parameters used in our calculations 
the primary cracks have an angular distribution with an 
average separation from 45 to 60 degrees. Therefore the 
mass of a fragment resulting from these plane meridional 
cracks is of the order of 10% of the sample mass, although 
typically only two to four cracks actually reach the sur- 
face breaking the material. This estimate corresponds 
to the range of masses that present the broad peak in 
the fragment mass distribution. This feature of the mass 
distribution function is not observed in the results of 2D 
simulations [HI |21] or 3D simulations of shell fragmenta- 
tion |491 152] , where obviously meridional cracks can not 
exist. 
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FIG. 10: (a) Fragment mass distribution for diflerent initial 
velocities. The straight line corresponds to a power-law with 
exponent -1.9 (b) Fragment size distribution weighted by mass 
for initial velocities with identical legend as above. 



VI. CONCLUSIONS 

We studied a brittle, disordered fragmenting solid 
sphere. We performed 3D DEM simulations with 3D 
beam-truss elements for the particle cohesion. Due to 
this computationally more laborious approach as com- 
pared to previous works, we were able to obtain a clearer 
picture of the fragmentation process, the evolution of 
fragmentation mechanisms, and its consequences for the 
fragment mass distribution. To get a clearer insight into 
the fracture initiation, we used continuum solutions for 
the stress field, obtained by the Finite Element Method. 
We were able to show, that 2D simulations for frag- 
menting systems are not capable of capturing fragmenta- 
tion by meridional cracks, that are the primary cracking 
mechanism. We found that cracks form inside the sample 
in the region above a compressive cone long time before 



FIG. 11: Fragment mass distribution at v = 122.5 m/s for 
different disorder in the bond breaking thresholds. The solid 
lines correspond to a power law with an exponential cuttoff 
for lower masses and the WeibuU distribution for large masses 
(Eq. 0) 



they are experimentally observable from the outside, if at 
all. They grow to form meridional fracture planes that 
result in a small number of large wedge shaped fragments, 
typically two to four. The increasing tensile radial and 
circumferential stress in the ring-shaped region above the 
contact plane gives rise to meridional cracks. The de- 
crease in the angular separation between these cracks 
could be explained by the Mott fragmentation model. 
Some of these cracks grow to form the meridional frac- 
ture planes that break the material in a small number of 
large fragments, and it is only then, that they become 
visible in experiments. 

The resulting mass distribution of the fragments 
presents a power law regime for small fragments and a 
broad peak for large fragments that can be fitted with 
a two-parameter Weibull distribution, in agreement with 
experimental results [TOl [I3J UHl [3D] • The fragment mass 
distribution is quite robust, independent on the macro- 
scopic material properties such as material strength and 
disorder distribution. Only the large fragment range of 
the mass distribution happens to be energy dependent, 
due to additional fragmentation processes that arise as 
one increases the impact energy. 

Even though our results are valid for various materi- 
als with disorder, they are limited to the class of brittle, 
heterogeneous media. Extensions to ductile materials are 
in progress. Another class of interesting questions deal 
with the problem of size effects, the influence of polydis- 
perse particles or the stiffness contrast of particles and 
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bcam-clcmcnts. The ability of the model to reproduce 
well defined crack planes also opens up the possibility 
to study other crack propagation problems in 3D. For 
technological applications questions about the influence 
of target shapes and the optimization potential to obtain 
desired fragment size distributions or to reduce impact 
energies are of broad interest. 
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